#include <iostream>
#include <fstream>

#include "bayesianclustering.h"

using namespace std;

//Global Variables
//*******************************
int N;
int Step;
int k=1;
double* Distrk;
double* Time;
double* Pi;
double* PiPrev;
int* Center;
int* Abundance;
double* Sigma;
bool* SigmaChange;
double* Likelihood;
double* LikelihoodRatio; //for UpdateLikelihood;
//************************************
double** LikelihoodTemp;
//************************************
char** SeqID;     //for Sequence ID
char** seq;  //for reading Fasta
char*** EigenSeqs=NULL; //for further iteration use
float** X;         //Distance Matrix
//double X2[MAXSize][MAXSize];
int* Z;       //missing data
int* NN;       //NN[i-1]=# of samples in i-th cluster
int* IterNN;
int* TempNN;
double* Centerlikelihood;  //square sum of current center
gsl_rng *r;
bool* Weight;  //for Birth-Death control;
unsigned long int WeightSum;
//*******************************

//Result Process
//*******************************
int Resultk;
double* ResultPi;
int* ResultCenter;
double* ResultSigma;
double* ResultLikelihood;
//*******************************

//Hyperparameters
//*******************************
//For Pi
double* Gamma;
//For Sigma
const double Alpha=2;
double Beta=3;
const double g=0.2;
double h;
double SigmaSum;
double Lower, Upper;
//For New-born Pi
double Gamma1=1;
//For Birth-rate
double Lamdab=1.0;
//For Stationary Time
double Lamda=2;
//Time Inspection
//*******************************
double TimeZ=0;
double TimeCenter=0;
double TimeSigma=0;
double TimePi=0;
double TimeBirth=0;
double TimeDeath=0;
double TimeLikelihood=0;

BayesianClusteringResult BayesianClustering(const char* fname, int step, bool ctrl, BayesianClusteringResult &CResult, int iter)
{   
    int len;
   	time_t start,stop;
    Step=step;
    k=1;
    if(ctrl==true)
    {
        Lower=1.5; Upper=2.5;
    }
    else
    {
        Lower=1; Upper=1.5;
    }
    InitRandomGenerator();
    LoadDistanceMatrix(fname,iter);
    GenerateInitial(iter);    
    //*******for time counting***********
	start=time(NULL);
    memset(Time,0,sizeof(double)*Step);
    memset(Distrk,0,sizeof(double)*Step);
    //***********************************
    cout<<Step<<endl;                
    cout<<"******Progress******"<<endl;
    cout<<"*";
    for(int i=0;i<Step;i++)
    {
            ProcessSampling(i,iter);
            if((20*i/Step)!=(20*(i-1)/Step))
                cout<<"*"; 
            /*             
            if(i>=2000)
            {
                //for result record
                CriteriaCalc(k);
                Pu=(double)(Pu*i+Purity)/(i+1);
                V=(double)(V*i+VI)/(i+1);
                ED=(double)(ED*i+EditDistance)/(i+1);           
            }
            */
    }
    cout<<endl;
    //for time counting
	stop=time(NULL);
	cout<<"Time Cost: "<<stop-start<<" sec"<<endl;
    //fs<<Pu<<","<<V<<","<<ED<<endl;
    cout<<Resultk<<endl;
    cout<<k<<endl;
    ResultProcess(); 
    CResult.ClusterNumber=k;
    int TempK=k;
    bool TempKFlag=false;
    if(k==1)
    {
        TempK++;
        TempKFlag=true;
    }
    CResult.ClusterCenterID=new char*[TempK];
    CResult.ClusterCenterSeq=new char*[TempK];
    CResult.ClusterStandardDeviation=new double[TempK];
    CResult.ClusterSize=new int[TempK];
    CResult.ClusterIterationSize=new int[TempK];
    CResult.ClusterMember=new string*[TempK];
    CResult.ClusterMemberSeq=new string*[TempK];
    CResult.ClusterMemberFlag=new bool*[TempK];
    TempNN=new int[TempK];
    if(TempKFlag==true)    //When Resultk=1
    {
        CResult.ClusterCenterID[TempK-1]=NULL;
        CResult.ClusterCenterSeq[TempK-1]=NULL;
        CResult.ClusterStandardDeviation[TempK-1]=0;
        CResult.ClusterSize[TempK-1]=0;
        CResult.ClusterIterationSize[TempK-1]=0;
        CResult.ClusterMember[TempK-1]=NULL;
        CResult.ClusterMemberSeq[TempK-1]=NULL;
        CResult.ClusterMemberFlag[TempK-1]=NULL;
        TempNN[TempK-1]=0;
    }
    for(int i=0;i<k;i++)
    {
        len=strlen(SeqID[Center[i]])+1;
        CResult.ClusterCenterID[i]=new char[len];
        memcpy(CResult.ClusterCenterID[i], SeqID[Center[i]], sizeof(char)*len);
        len=strlen(seq[Center[i]])+1;
        CResult.ClusterCenterSeq[i]=new char[len];
        memcpy(CResult.ClusterCenterSeq[i], seq[Center[i]], sizeof(char)*len);
        CResult.ClusterStandardDeviation[i]=Sigma[i];
        CResult.ClusterSize[i]=NN[i];
        CResult.ClusterIterationSize[i]=IterNN[i];
        CResult.ClusterMember[i]=new string[IterNN[i]];
        CResult.ClusterMemberFlag[i]=new bool[IterNN[i]];
        memset(CResult.ClusterMemberFlag[i],0,sizeof(bool)*IterNN[i]);
        CResult.ClusterMemberSeq[i]=new string[MIN(20,NN[i])];
    }
    memset(TempNN, 0,  sizeof(int)*k);   
    for(int i=0;i<N;i++)
    {
        CResult.ClusterMember[Z[i]][TempNN[Z[i]]].assign(SeqID[i]);
        TempNN[Z[i]]++;
    }
    for(int i=0;i<k;i++)
    {
        if(TempNN[i]!=IterNN[i])
            cout<<"Unexpected Result Error in Block"<<endl;
    }
    cout<<"Basic Operations Done"<<endl;
    //*****************************************Generate Eigen Seqs*********************************************************
    memset(TempNN,0,sizeof(int)*k);
    int TempSeqNum, temp;
    for(int i=0;i<N;i++)
    {
        TempNN[Z[i]]+=MIN(20,Abundance[i]);
    }
    if(iter!=0)
    {
    for(int i=0;i<k;i++)
    {

        unsigned long int Range=MIN(20,TempNN[i]);
        if(Range==20)
        {
            for(int j=0;j<20;j++)
            {
                TempSeqNum=gsl_rng_uniform_int(r, TempNN[i])+1;
                temp=0;
                for(int l=0;l<N;l++)
                {
                    if(Z[l]==i)
                    {
                        temp+=MIN(20,Abundance[l]);
                        if(temp>=TempSeqNum)
                        {
                            CResult.ClusterMemberSeq[i][j].assign(EigenSeqs[l][TempSeqNum-1-(temp-MIN(20,Abundance[l]))]);
                            break;
                        }
                    }
                }
            }
        }
        else
        {
            int TempCount=0;
            for(int j=0;j<N;j++)
            {
                if(Z[j]==i)
                {
                    for(int l=0;l<MIN(20,Abundance[j]);l++)
                    {
                        CResult.ClusterMemberSeq[i][TempCount].assign(EigenSeqs[j][l]);
                        TempCount++;
                    }
                }
            }
        }
    }
    }//if iter!=0
    
    else
    {
        for(int i=0;i<k;i++)
        {
            for(int j=0;j<MIN(20,TempNN[i]);j++)
            {
               TempSeqNum=gsl_rng_uniform_int(r, TempNN[i]);
                temp=0;
                for(int l=0;l<N;l++)
                {
                    if(Z[l]==i)
                    {
                        temp+=MIN(20,Abundance[l]);
                        if(temp>=TempSeqNum)
                        {
                            CResult.ClusterMemberSeq[i][j].assign(seq[l]);
                            break;
                        }//if temp>TempSeqNum
                    }//if Z[l]==i
                }//l
            }//j
        }//i
    }//if iter==0
    cout<<"Choosing Seqs Compelete"<<endl;
    //********************************************************************************************************************
            
    DELETE(TempNN);
    //cout<<TimeZ<<" "<<TimeCenter<<" "<<TimePi<<" "<<TimeSigma<<" "<<TimeBirth<<" "<<TimeDeath<<" "<<TimeLikelihood<<endl;
    //cout<<"Build Complete"<<endl;
    gsl_rng_free(r);
    DeleteMem(N,iter);
    return CResult;    
}

void InitRandomGenerator()
{
     r = gsl_rng_alloc (gsl_rng_taus);
     gsl_rng_set(r,(unsigned int)time((time_t *)NULL));
     return;   
}

void LoadDistanceMatrix(const char* fname, int iter)
{
    fstream tempout1,tempout2;
    float maxtemp=0.0, mintemp=100.0;
    int matrix = EDNAFULL;
	unsigned int seq_id = 0, tar_id = 0;
	alignment (*alignTool) (const char *, const char *, alignment, int, int) = NULL, align = {0, 0, 0};
	unsigned int UnN=(unsigned int)N;
	char* ffname;
	ffname=new char[strlen(fname)];
    strcpy(ffname, fname);    
    if(iter==0)
    {
    seq = readFastASeq(ffname, &UnN);
    N=(int)UnN; 
    NewMem(N);
    for(int i=0;i<N;i++)
        Abundance[i]=1; 
    alignTool = NW_alignmentEx;
	for(seq_id = 0; seq_id < N-1; seq_id++){
        X[seq_id][seq_id]=0.0;
		for(tar_id = seq_id+1; tar_id < N; tar_id++){
			align = alignTool(seq[seq_id], seq[tar_id], align, matrix, 1);
			X[seq_id][tar_id]=(float)(1.0-(double)align.matches/(double)SHORT_ALIGN(align))*100.0;
		}
	}
	X[N-1][N-1]=0.0;
    for(seq_id=0; seq_id < N; seq_id++)
        for(tar_id = 0; tar_id < seq_id; tar_id++)
            X[seq_id][tar_id]=X[tar_id][seq_id];
    }//if iter==0
    
    else
    {   
    seq = readBlockSeq(ffname, &UnN);
    N=(int)UnN;
    //k=N;   
    alignTool = NW_alignmentEx;
    if(iter==20)
    {
        tempout1.open("tempmatrix1",ios_base::out);
        tempout2.open("tempmatrix2",ios_base::out);
    }
	for(seq_id = 0; seq_id < N; seq_id++){
        X[seq_id][seq_id]=Sigma[seq_id];
        if(iter==20)
        {
            tempout1<<SeqID[seq_id]<<" ";
            tempout2<<SeqID[seq_id]<<" ";
        }
	    for(tar_id = 0; tar_id < N; tar_id++){
            if(tar_id>seq_id)
            {
			align = alignTool(seq[seq_id], seq[tar_id], align, matrix, 1);
			X[seq_id][tar_id]=(float)(1.0-(double)align.matches/(double)SHORT_ALIGN(align))*100.0;
            X[tar_id][seq_id]=X[seq_id][tar_id];
            if(iter==20)
            {
                tempout1<<X[seq_id][tar_id]<<" ";
                if(X[seq_id][tar_id]>=15)
                    tempout2<<X[seq_id][tar_id]<<" ";
            }
			if((X[seq_id][tar_id]<15)&&(X[seq_id][tar_id]>(float)Lower/2))
			{
                for(int i=0;i<MIN(20,Abundance[seq_id]);i++)
                {
                    align = alignTool(seq[tar_id], EigenSeqs[seq_id][i], align, matrix, 1);
                    X[seq_id][tar_id]+=(float)(1.0-(double)align.matches/(double)SHORT_ALIGN(align))*100.0;
                }
                X[seq_id][tar_id]/=(MIN(20,Abundance[seq_id])+1);
                for(int i=0;i<MIN(20,Abundance[tar_id]);i++)
                {
                    align = alignTool(seq[seq_id], EigenSeqs[tar_id][i], align, matrix, 1);
                    X[tar_id][seq_id]+=(float)(1.0-(double)align.matches/(double)SHORT_ALIGN(align))*100.0;
                }
                X[tar_id][seq_id]/=(MIN(20,Abundance[tar_id])+1);
                if(iter==20)
                    tempout2<<X[seq_id][tar_id]<<" ";
            }
            }//if
            if((tar_id<=seq_id)&&(iter==20))
            {
                tempout1<<X[seq_id][tar_id]<<" ";
                tempout2<<X[tar_id][seq_id]<<" ";
            }
		}
		if(iter==20)
	    {
                tempout1<<endl;
                tempout2<<endl;            
        }
        /*
		for(tar_id = seq_id+1; tar_id < N; tar_id++){
			align = alignTool(seq[seq_id], seq[tar_id], align, matrix, 1);
			X[seq_id][tar_id]=(float)(1.0-(double)align.matches/(double)SHORT_ALIGN(align))*100.0;
            X[tar_id][seq_id]=X[seq_id][tar_id];

			if((X[seq_id][tar_id]<15)&&(X[seq_id][tar_id]>Upper/2))
			{
                for(int i=0;i<MIN(20,Abundance[seq_id]);i++)
                {
                    align = alignTool(seq[tar_id], EigenSeqs[seq_id][i], align, matrix, 1);
                    X[seq_id][tar_id]+=(float)(1.0-(double)align.matches/(double)SHORT_ALIGN(align))*100.0;
                }
                X[seq_id][tar_id]/=(MIN(20,Abundance[seq_id])+1);
                for(int i=0;i<MIN(20,Abundance[tar_id]);i++)
                {
                    align = alignTool(seq[seq_id], EigenSeqs[tar_id][i], align, matrix, 1);
                    X[tar_id][seq_id]+=(float)(1.0-(double)align.matches/(double)SHORT_ALIGN(align))*100.0;
                }
                X[tar_id][seq_id]/=(MIN(20,Abundance[tar_id])+1);
            }
		}
		*/
	}
	X[N-1][N-1]=Sigma[N-1];
    }
    for(int i=0;i<N;i++)
    {   
        for(int j=0;j<N;j++)
        {
            if(X[i][j]>=maxtemp)
                 maxtemp=X[i][j];
            if((X[i][j]<=mintemp)&&(X[i][j]!=0))
                 mintemp=X[i][j];
        }
    }
    h=(double)100*g/(maxtemp-mintemp)/(maxtemp-mintemp)/Alpha;
    if(iter==20)
    {
        tempout1.close();
        tempout2.close();
    }
    delete[] ffname;
    return;    
}

int GenerateInitial(int iter)
{
    //if(iter==0)
    //{    
    //Init Weight & WeightSum    
    for(int i=0;i<N;i++)
            Weight[i]=true; //false
    memset(PiPrev,0,sizeof(double)*N);
    memset(Pi,0,sizeof(double)*N);
    WeightSum=N; //0
    //Init Pi
    for(int i=0;i<N;i++)
            Gamma[i]=1.0;
            //Pi[i]=0.2;
    gsl_ran_dirichlet(r, k, Gamma, Pi);
    memcpy(PiPrev, Pi, sizeof(double)*k);
    memcpy(ResultPi, Pi, sizeof(double)*k);
    //Init Center
    for(int i=0;i<k;i++)
    {
            //Center[i]=(int)gsl_ran_flat(r,100*i,100*i+99);
            Center[i]=(int)gsl_ran_flat(r,0,N);
            ResultCenter[i]=Center[i];
            Weight[Center[i]]=false;
            WeightSum--;
    }
    //Init Sigma
    SigmaSum=0;
    Beta=gsl_ran_gamma(r,g,(double)1/h);
    for(int i=0;i<k;i++)
    {
            for(int j=0;j<10;j++)
            {
                Sigma[i]=sqrt(1/gsl_ran_gamma(r, Alpha, (double)1/Beta));
                if((Sigma[i]<=Upper)&&(Sigma[i]>=Lower))
                    break;
            }
            while((Sigma[i]>Upper)||(Sigma[i]<Lower))
            {
                Sigma[i]=sqrt(1/gsl_ran_gamma(r, Alpha+1, (double)1/(Upper+Lower)));
            }
            ResultSigma[i]=Sigma[i];
            SigmaSum+=(double)1/Sigma[i]/Sigma[i];
    }
    memset(SigmaChange,0,sizeof(bool)*k);
    //Init LikelihoodRatio
    UpdateLikelihoodRatio();
    for(int i=0;i<N;i++)
        ResultLikelihood[i]=Likelihood[i];
    Resultk=k;
    //}//if(iter==0)
    /*
    else
    {
    //Init Weight & WeightSum    
    for(int i=0;i<N;i++)
            Weight[i]=false; //false
    WeightSum=0; //0
    //Init Pi
    for(int i=0;i<N;i++)
            Gamma[i]=1.0;
    double* GammaTemp=new double[k+1];
    memset(GammaTemp,0,sizeof(double)*(k+1));
    for(int i=0;i<k;i++)
        GammaTemp[i]=Gamma[i]+Abundance[i];
    memset(PiPrev,0,sizeof(double)*N);
    memset(Pi,0,sizeof(double)*N);
    gsl_ran_dirichlet(r, k, GammaTemp, Pi);      
    delete[] GammaTemp;
    memcpy(PiPrev, Pi, sizeof(double)*k);
    memcpy(ResultPi, Pi, sizeof(double)*k);
    //Init Center
    for(int i=0;i<k;i++)
    {
            Center[i]=i;
            ResultCenter[i]=Center[i];
    }
    //Init Sigma
    SigmaSum=0;
    Beta=gsl_ran_gamma(r,g,(double)1/h);
    for(int i=0;i<k;i++)
    {
            ResultSigma[i]=Sigma[i];
            SigmaSum+=(double)1/Sigma[i]/Sigma[i];
    }
    memset(SigmaChange,0,sizeof(bool)*k);
    //Init LikelihoodRatio
    UpdateLikelihoodRatio();
    for(int i=0;i<N;i++)
        ResultLikelihood[i]=Likelihood[i];
    delete[] GammaTemp;
    Resultk=k;        
    }//else
    */
    return 1;   
}

int UpdateZ()
{
     time_t start;
     start=time(NULL);
     memset(NN,0,sizeof(int)*k);
     double z,temp=0;
     bool flag;
     //******************************************
     for(int i=0;i<N;i++)
     {
          flag=false;
          for(int j=0;j<k;j++)
          {
                  if(i==Center[j])
                  {
                       Z[i]=j;
                       NN[j]+=Abundance[i];
                       flag=true;
                       break;
                  }
                       
          }
          if(flag==false)
          {   
             z=gsl_ran_flat(r, 0, Likelihood[i]);
             temp=0;
             for(int j=0;j<k;j++)
             {
                     temp+=LikelihoodTemp[i][j];
                     if(temp>=z)
                     {
                         Z[i]=j;
                         NN[j]+=Abundance[i];
                         break;
                     }
             }
          }
     }
     /*
     if(k==1)
     {
         for(int i=0;i<N;i++)
             if(Z[i]!=0)
                 cout<<"Z["<<i<<"]";
     }
     */
     //*******************************************
     time_t stop;
     stop=time(NULL);
     TimeZ+=stop-start;
     //cout<<"Z";
     return 1;
}

int UpdatePi()
{
     time_t start;
     start=time(NULL);
       int d=0;
       //Sampling Method
       memcpy(PiPrev,Pi,sizeof(double)*k);
       double* GammaTemp=new double[k+1];
       memset(GammaTemp,0,sizeof(double)*(k+1));
       for(int i=0;i<k;i++)
               GammaTemp[i]=Gamma[i]+NN[i];
       gsl_ran_dirichlet(r, k, GammaTemp, Pi);      
       //Posterior Mean Method        
       /*
       for(int i=0;i<k;i++)
               d+=Gamma[i]+NN[i];
       for(int i=0;i<k;i++)
               Pi[i]=(double)(Gamma[i]+NN[i])/d;
       */
       delete[] GammaTemp;
     time_t stop;
     stop=time(NULL);
     TimePi+=stop-start;  
     //cout<<"P";            
       return 1;
}

int UpdateCenter(int m, bool f)
{
     time_t start;
     start=time(NULL);
     double *p=NULL,*rand=NULL,*temp=NULL;
     bool *flag=NULL;
     if(m%10==0)
     {
     memset(Centerlikelihood,0,sizeof(double)*k);
     for(int i=0;i<N;i++)
     {
         double x=X[i][Center[Z[i]]];
         Centerlikelihood[Z[i]]+=((x*x/2)*Abundance[i]);
     }
     double* TempLikelihood=new double[N];      //for UpdateCenter;
     memset(TempLikelihood, 0,sizeof(double)*N);
     //cout<<"C1";
     if(k>1)
     {
         p=new double[k];
         rand=new double[k];
         temp=new double[k];
         flag=new bool[k];
         memset(rand,0,sizeof(double)*k);
         memset(flag,0,sizeof(bool)*k);
         memset(p,0,sizeof(double)*k);
         memset(temp,0,sizeof(double)*k);
     }
     else
     {
         p=new double[k+1];
         rand=new double[k+1];
         temp=new double[k+1];
         flag=new bool[k+1];
         memset(rand,0,sizeof(double)*(k+1));
         memset(flag,0,sizeof(bool)*(k+1));
         memset(p,0,sizeof(double)*(k+1));
         memset(temp,0,sizeof(double)*(k+1));
     }
     //cout<<"C2";
     for(int i=0;i<N;i++)
     {
             for(int j=0;j<N;j++)
             {
                  if(Z[j]==Z[i])
                  {
                      double x=X[j][i];
                      TempLikelihood[i]+=((x*x/2)*Abundance[j]);
                  }
             }
             TempLikelihood[i]-=Centerlikelihood[Z[i]];
             TempLikelihood[i]=gsl_expm1(-TempLikelihood[i]/Sigma[Z[i]]/Sigma[Z[i]])+1;
             p[Z[i]]+=MIN(10e300,TempLikelihood[i]);
     }
     //cout<<"C6";
     for(int j=0;j<k;j++)
             rand[j]=gsl_ran_flat(r, 0, p[j]);
     //cout<<"C7";
     for(int i=0;i<N;i++)
     {
             temp[Z[i]]+=MIN(10e300,TempLikelihood[i]);
             if((temp[Z[i]]>=rand[Z[i]])&&(flag[Z[i]]==0))
             {
                     Weight[Center[Z[i]]]=true;                                     
                     Center[Z[i]]=i;
                     Weight[i]=false;
                     flag[Z[i]]=1;
             }
     }
     //cout<<"C3";
    delete[] temp;
    delete[] rand;
    delete[] p;
    delete[] TempLikelihood;
    delete[] flag;
   // cout<<"C4";
    time_t stop;
    stop=time(NULL);
    TimeCenter+=stop-start;
    return 1;
    }
    else if(f==1)
    {
     Centerlikelihood[k-1]=0;
     for(int i=0;i<N;i++)
         if(Z[i]==k-1)
         {
             double x=X[i][Center[k-1]];
             Centerlikelihood[k-1]+=((x*x/2)*Abundance[i]);
         }
     double* TempLikelihood=new double[N];      //for UpdateCenter;
     memset(TempLikelihood, 0,sizeof(double)*N);
     if(k>1)
     {
         p=new double[k];
         rand=new double[k];
         temp=new double[k];
         flag=new bool[k];
         memset(rand,0,sizeof(double)*k);
         memset(flag,0,sizeof(bool)*k);
         memset(p,0,sizeof(double)*k);
         memset(temp,0,sizeof(double)*k);
     }
     else
     {
         p=new double[k+1];
         rand=new double[k+1];
         temp=new double[k+1];
         flag=new bool[k+1];
         memset(rand,0,sizeof(double)*(k+1));
         memset(flag,0,sizeof(bool)*(k+1));
         memset(p,0,sizeof(double)*(k+1));
         memset(temp,0,sizeof(double)*(k+1));
     }
     for(int i=0;i<N;i++)
     {
         if(Z[i]==k-1)
         {
             for(int j=0;j<N;j++)
             {
                  if(Z[j]==k-1)
                  {
                      double x=X[j][i];
                      TempLikelihood[i]+=((x*x/2)*Abundance[j]);
                  }
             }
             TempLikelihood[i]-=Centerlikelihood[k-1];
             TempLikelihood[i]=gsl_expm1(-TempLikelihood[i]/Sigma[k-1]/Sigma[k-1])+1;
             p[k-1]+=TempLikelihood[i];
         }
     }
     rand[k-1]=gsl_ran_flat(r, 0, p[k-1]);
     for(int i=0;i<N;i++)
     {
         if(Z[i]==k-1)
         {
             temp[k-1]+=TempLikelihood[i];
             if((temp[k-1]>=rand[k-1])&&(flag[k-1]==0))
             {
                     Weight[Center[k-1]]=true;                                     
                     Center[k-1]=i;
                     Weight[i]=false;
                     flag[k-1]=1;
             }
        }
    }
   // cout<<"C3";
    delete[] temp;
    delete[] rand;
    delete[] p;
    delete[] TempLikelihood;
    delete[] flag;
    }
    //cout<<"C4";
    time_t stop;
    stop=time(NULL);
    TimeCenter+=stop-start;
    return 1;      
}

int UpdateSigma()
{
     time_t start;
     start=time(NULL);
     double temp;
    Beta=gsl_ran_gamma(r, (double)g+k*Alpha, (double)1/(h+SigmaSum));
    SigmaSum=0;
    memset(SigmaChange,0,sizeof(bool)*k);
    //double temp=0;
    //for(int i=0;i<k;i++)
        //temp+=Centerlikelihood[i]; 
    for(int i=0;i<k;i++)
    {
            //Sigma[i]=sqrt(1/gsl_ran_gamma(r, Alpha+(double)NN[i]/2, (double)1/(Beta+temp)));   
            temp=sqrt(1/gsl_ran_gamma(r, Alpha+(double)NN[i]/2, (double)1/(Beta+Centerlikelihood[i])));
            /*   
            if(Sigma[i]>Threshold)
                Sigma[i]=sqrt(1/gsl_ran_gamma(r, Alpha*10+1, (double)1/Threshold/10));
            if(Sigma[i]<1)
                Sigma[i]=sqrt(1/gsl_ran_gamma(r, Alpha*10+1, (double)1/Threshold/10));
            */
            /*
            while((Sigma[i]>Threshold)||(Sigma[i]<1))
                Sigma[i]=sqrt(1/gsl_ran_gamma(r, Alpha+1, (double)1/Threshold));
            */
           while((temp>Upper)||(temp<Lower))
                temp=sqrt(1/gsl_ran_gamma(r, Alpha+1, (double)1/(Upper+Lower)));
           double ratio=temp/Sigma[i];
           if((ratio<1.2)&&(ratio>0.8))
                SigmaChange[i]=true;
           else
               Sigma[i]=temp;
           SigmaSum+=1/Sigma[i]/Sigma[i];    
    }
     time_t stop;
     stop=time(NULL);
     TimeSigma+=stop-start;
    // cout<<"S";
    return 1;
}

int Birth(int m)
{
     time_t start;
     start=time(NULL);
    //Init new state
    Pi[k-1]=gsl_ran_beta(r, Gamma1, (double)k*Gamma1);
    //**********************choose intial sigma*****************************
    for(int j=0;j<10;j++)
    {
        Sigma[k-1]=sqrt(1/gsl_ran_gamma(r, Alpha, (double)1/Beta));
        if((Sigma[k-1]<=Upper)&&(Sigma[k-1]>=Lower))
            break;
    }
    while((Sigma[k-1]>Upper)||(Sigma[k-1]<Lower))
    {
        Sigma[k-1]=sqrt(1/gsl_ran_gamma(r, Alpha+1, (double)2/(Upper+Lower)));
    }
    //************************************************************************ 
    for(int i=0;i<k-1;i++)
    {
        Pi[i]*=(1-Pi[k-1]);
        for(int j=0;j<N;j++)
            LikelihoodTemp[j][i]*=(1-Pi[k-1]);
    } 
    for(int i=0;i<N;i++)
    {
        LikelihoodTemp[i][k-1]=Pi[k-1]*gaussian(X[i][Center[k-1]],Sigma[k-1]);
        Likelihood[i]=Likelihood[i]*(1-Pi[k-1])+LikelihoodTemp[i][k-1];
    }
    //Update whole parameters
    //cout<<"Birth "<<Center[k-1]<<endl;
     time_t stop;
     stop=time(NULL);
     TimeBirth+=stop-start;
    // cout<<"B";
    UpdateZ();
    UpdatePi();
    UpdateCenter(m,1);
    UpdateSigma();
    UpdateLikelihoodRatio();    
    return 1;       
}

int Death(int j,int m)
{
     time_t start;
     start=time(NULL);
    //Init new state
    for(int i=0;i<N;i++)
            Likelihood[i]=(Likelihood[i]-LikelihoodTemp[i][j])/(1-Pi[j]);
    for(int i=0;i<k;i++)
    {
            if(i!=j)
                Pi[i]/=(1-Pi[j]);
            for(int l=0;l<N;l++)
            {
                if(i<j)
                    LikelihoodTemp[l][i]/=(1-Pi[j]);
                else
                   if(i<k-1)
                    LikelihoodTemp[l][i]=LikelihoodTemp[l][i+1]/(1-Pi[j]);
            }
    }
    Weight[Center[j]]=true;
    WeightSum++;
    k--;
    for(int i=j;i<k;i++)
    {
            Pi[i]=Pi[i+1];
            Center[i]=Center[i+1];
            Sigma[i]=Sigma[i+1];
    }
     time_t stop;
     stop=time(NULL);
     TimeDeath+=stop-start;
    //Update whole Parameters
    //cout<<"Death "<<j<<endl;
   // cout<<"D";
    UpdateZ();
    UpdatePi();
    UpdateCenter(m,0);
    UpdateSigma(); 
    UpdateLikelihoodRatio();    
    return 1;    
}

int ProcessSampling(int m, int iter)
{
    bool flag=false;
    double temp=0.0,z,Ti,L=0.0;
    flag=false;
    if(k==1)
    {
           if(m==0)
               Ti=0;
           else
               Ti=Time[m-1];
           Time[m]=gsl_ran_exponential(r,(double)1/Lamdab);
           Distrk[k]+=Time[m];
           for(int i=0;i<N;i++)
               if(Weight[i]==true)
               {         
                   //if(Likelihood[i]>=10.0e299||Likelihood[i]<=10.0e-300)
                       //cout<<Likelihood[i]<<" ";
                   temp+=1.0/MAX(10e-300,Likelihood[i]);
                   //if(Likelihood[i]>=10.0e299||Likelihood[i]<=10.0e-300)
                       //cout<<temp<<" ";
               }
           //cout<<temp;      
           z=gsl_ran_flat(r, 0, temp);
           temp=0;
           for(int i=0;i<N;i++)
           {
               if(Weight[i]==true)
               {
                   temp+=1.0/MAX(10.0e-300,Likelihood[i]); 
                   if(temp>=z)
                   {
                       Center[k]=i;
                       Weight[i]=false;
                       WeightSum--;
                       k++;
                       flag=true;
                       break;
                   }
               }
           }
           //cout<<"[1,"<<k<<"]";
           Birth(m);
           return 1;
    }
    else
    {
        if(m==0)
            Ti=0;
        else
           Ti=Time[m-1];
        for(int i=0;i<k;i++)
                temp+=LikelihoodRatio[i];
        if(k<N)
                temp+=Lamdab;//(double)k;
        Time[m]=gsl_ran_exponential(r,(double)1/temp);
        Distrk[k]+=Time[m];
        z=gsl_ran_flat(r, 0, temp);
        temp=0;
        for(int i=0;i<k;i++)
        {
                temp+=LikelihoodRatio[i];
                if(temp>=z)
                {
                           //cout<<"2";
                           Death(i,m);
                           flag=true;
                           break;
                }
        }
        if(flag==false)
        {
           temp=0;
       	       //Randomly choose a new group center from data.
           for(int i=0;i<N;i++)
               if(Weight[i]==true)
                   temp+=1.0/MAX(10.0e-300,Likelihood[i]);        
           z=gsl_ran_flat(r, 0, temp);
           temp=0;
           for(int i=0;i<N;i++)
           {
               if(Weight[i]==true)
               {
                       temp+=1.0/MAX(10.0e-300,Likelihood[i]);
                   if(temp>=z)
                   {
                       Center[k]=i;
                       Weight[i]=false;
                       WeightSum--;
                       k++;
                       flag=true;
                       break;
                   }
               }
           }
           if(flag==true)
           {
              //cout<<"[3,"<<k<<"]";
              Birth(m);
           }
        }
    }
    for(int i=0;i<N;i++)
    {
        if(Likelihood[i]!=0&&ResultLikelihood[i]==0)
            L+=log10(Likelihood[i])+300;
        if(Likelihood[i]==0&&ResultLikelihood[i]==0)
            continue;
        if(Likelihood[i]!=0&&ResultLikelihood[i]!=0)
            L+=log10(Likelihood[i])-log10(ResultLikelihood[i]);
        if(Likelihood[i]==0&&ResultLikelihood[i]!=0)
            L-=300;
        /*
        if(Likelihood[i]==0)
            if(ResultLikelihood[i]!=0)
                L*=(10E-308/ResultLikelihood[i]);
                else
                    L*=1.0;
            else
                if(ResultLikelihood[i]!=0)
                    L*=(Likelihood[i]/ResultLikelihood[i]);
                else
                    L*=(Likelihood[i]/10E-308);
        */
    }
    if(L>=307)
        temp=10.0E307;
    if(L<=-308)
        temp=0;
    if(L<307&&L>-308)
        temp=pow((double)10,L);
    L=temp;
        //if((L>1)&&(temp<MINProb+1))
    if(L>pow(Lamdab,Resultk-k)*factorial(Resultk,k))
    {
        if(iter==20)
        {
            fstream TempRatioMoniter;
            TempRatioMoniter.open("LikelihoodRatio.txt",ios_base::out|ios_base::app);
            if(TempRatioMoniter.fail())
                cout<<"Cannot write likelihoodratio info!"<<endl;
            for(int i=0;i<k;i++)
            {
                TempRatioMoniter<<LikelihoodRatio[i]<<" ";
                //cout<<LikelihoodRatio[i]<<" ";
            }
            TempRatioMoniter<<endl;
            //cout<<endl;
            TempRatioMoniter.close();
        }                                            
        Resultk=k;
        for(int i=0;i<k;i++)
        {
            ResultPi[i]=Pi[i];
            ResultCenter[i]=Center[i];
            ResultSigma[i]=Sigma[i];
        }
        for(int i=0;i<N;i++)
            ResultLikelihood[i]=Likelihood[i];
    }
    return 1;       
}

int UpdateLikelihoodRatio()
{
    time_t start;
    start=time(NULL);
    memset(Likelihood, 0, sizeof(double)*N);
    for(int i=0;i<N;i++)
            for(int j=0;j<k;j++)
            {
                 if(SigmaChange[j]==false)
                     LikelihoodTemp[i][j]=Pi[j]*gaussian(X[i][Center[j]],Sigma[j]);
                 else
                     LikelihoodTemp[i][j]*=(Pi[j]/PiPrev[j]);                  
                 Likelihood[i]+=LikelihoodTemp[i][j]; 
            }
/*
    for(int i=0;i<N;i++)
            if(Likelihood[i]==0)
                Likelihood[i]=4e-300;
*/
    if(k==1)
        LikelihoodRatio[0]=0;
    else
    {
        double* TempLikelihood=new double[N];
        double* LogLikelihoodRatio=new double[k+1];
        memset(LogLikelihoodRatio,0,sizeof(double)*(k+1));     
        for(int j=0;j<k;j++)
        {
            LikelihoodRatio[j]=1.0;
            for(int i=0;i<N;i++)
            {
                    if(Likelihood[i]==0)
                         TempLikelihood[i]=1;
                    else
                         TempLikelihood[i]=(Likelihood[i]-LikelihoodTemp[i][j])/(1-Pi[j])/Likelihood[i];
                    LogLikelihoodRatio[j]+=Abundance[i]*log10(TempLikelihood[i]);
                    //LikelihoodRatio[j]*=pow(TempLikelihood[i],Abundance[i]);
                    //if(LikelihoodRatio[j]<0)
                        //cout<<pow(TempLikelihood[i],Abundance[i])<<endl;
            }
            if(LogLikelihoodRatio[j]<=-308)
                LikelihoodRatio[j]=0;
            if(LogLikelihoodRatio[j]>=307)
                LikelihoodRatio[j]=10.0e300;
            if((LogLikelihoodRatio[j]<307)&&(LogLikelihoodRatio[j]>-308))
                LikelihoodRatio[j]=pow((double)10,LogLikelihoodRatio[j]);
        }
        delete[] TempLikelihood;
        delete[] LogLikelihoodRatio;
    }
     time_t stop;
     stop=time(NULL);
     TimeLikelihood+=stop-start;
    // cout<<"L";
    return 1;                       
}

int ResultProcess()
{
    
     if(Resultk==1)
     {
         Pi[0]=ResultPi[0];
         Center[0]=ResultCenter[0];
         Sigma[0]=ResultSigma[0];
         IterNN[0]=0;
         NN[0]=0;
         k=Resultk;
         for(int i=0;i<N;i++)
         {
             Z[i]=0;
             NN[0]+=Abundance[i];
             IterNN[0]++;
         }
     }
     
     
     else
     {
     memcpy(Pi, ResultPi, sizeof(double)*Resultk);
     memcpy(Center, ResultCenter, sizeof(int)*Resultk);
     memcpy(Sigma, ResultSigma, sizeof(double)*Resultk);
     memset(NN,0,sizeof(int)*Resultk);
     memset(IterNN,0,sizeof(int)*Resultk);
     k=Resultk;
     double z,temp,maxtemp;
     int Zmax;
     bool flag;
     for(int i=0;i<N;i++)
     {
          maxtemp=0;
          
          flag=false;
          for(int j=0;j<k;j++)
          {
                  if(i==Center[j])
                  {
                       Z[i]=j;
                       NN[j]+=Abundance[i];
                       IterNN[j]++;
                       flag=true;
                       break;
                  }                       
          }
          if(flag==false)
          {   
             for(int j=0;j<k;j++)
             {                     
                     temp=Pi[j]*gaussian(X[i][Center[j]],Sigma[j]);
                     if(temp>maxtemp)
                     {
                         maxtemp=temp;
                         Zmax=j;
                     }
             }
             Z[i]=Zmax;
             NN[Zmax]+=Abundance[i];
             IterNN[Zmax]++;
          }
    }
    }//else
   // cout<<"R";        
    return 1;   
}

/*
int CriteriaCalc(int t)
{
    int temp,tmp;
    //Purity
    Purity=(double)t;
    int Candidate[MAXSize];
    memset(Candidate, -1, sizeof(Candidate));
    for(int i=0;i<t;i++)
        for(int j=0;j<N;j++)
        {
                if(Z[j]==i)
                {
                    if((Candidate[i]!=-1)&&(StandardZ[j]!=Candidate[i]))
                    {
                            Purity--;
                            break;
                    }
                    else
                         Candidate[i]=StandardZ[j];
                }
        }
    tmp=t;
    for(int i=0; i<t; i++)
        if(NN[i]==0)
        {
            Purity--;
            tmp--;
        }
    Purity/=(double)tmp;
    
    //VI
    double HResult=0.0;
    double HStandard=0.0;
    double HRandS=0.0;
    for(int i=0;i<t;i++)
    {
        if(NN[i]!=0)
            HResult-=(double)NN[i]/N*log((double)NN[i]/N);
    }
    for(int i=0;i<Standardk;i++)
        HStandard-=(double)StandardNN[i]/N*log((double)StandardNN[i]/N);
    for(int i=0;i<t;i++)
        for(int j=0;j<Standardk;j++)
        {
            temp=0;
            for(int l=0;l<N;l++)
            {
                    if((Z[l]==i)&&(StandardZ[l]==j))
                        temp++;
            }
            if(temp>0)
                HRandS+=(double)temp/N*log((double)temp/NN[i]/StandardNN[j]*N);
        }
    VI=HResult+HStandard-2*HRandS;
    //Edit Distance
    EditDistance=0;
    for(int i=0;i<t;i++)
        for(int j=0;j<Standardk;j++)
        {
            temp=0;
            for(int l=0;l<N;l++)
            {
                    if((Z[l]==i)&&(StandardZ[l]==j))
                    {
                        temp=1;
                        EditDistance++;
                        break;
                    }
            }
        }
    EditDistance*=2;
    EditDistance-=(tmp+Standardk);
    
    return 1;        
}
*/

double factorial(int t, int l)
{
    if(t<l)
        return l*factorial(t,l-1);
    if(t>l)
        return (double)1/t*factorial(t-1,l);
    if(t==l)
        return 1;
}

double gaussian(double x, double y)
{
    if(x>15)
        return 0;
    else
        return gsl_ran_gaussian_pdf(x,y);
}

char** readFastASeq(char *fname, unsigned int *num_seq)
{
	char dummy[MAX_LINE], **seq = NULL, tempdummy[MAX_LINE];
	unsigned int i = 0, len = 0, lentotal = 0;
	FILE *fp = fopen(fname, "r");

	if(fp == NULL) printf("Cannot open `%s' file.\n", fname);
	
    *num_seq = 0;
	while(feof(fp) == 0){
		fgets(dummy, MAX_LINE, fp);
		if(dummy[0] == '>')
			(*num_seq)++;
	}
	printf("num_seq = %u \n", *num_seq);
	rewind(fp);
    
	MALLOC(seq, char*, *num_seq);
	MALLOC(SeqID, char*, *num_seq);

	fgets(dummy, MAX_LINE, fp); /* Read first line */
	for(i = 0; i < *num_seq; i++){
		if(dummy[0] == '>'){
            len = strlen(dummy);
            len += (dummy[len-1] == '\n')? -1:0;
            MALLOC(SeqID[i], char, len+1);
            memcpy(SeqID[i], dummy, sizeof(char)*len);        
			if(fgets(dummy, MAX_LINE, fp) == NULL) printf("This is not a valid fasta file (no read for a header).\n");
			lentotal=0;   
			len = strlen(dummy);
			if(len == MAX_LINE) printf("There is very long line (> %u).\n", MAX_LINE);
			len += (dummy[len-1] == '\n')? -1:0;
            memcpy(tempdummy+lentotal, dummy, sizeof(char)*len);			
			lentotal+=len;
			if(fgets(dummy, MAX_LINE, fp) == NULL){ //end of file
			    MALLOC(seq[i], char, lentotal+1);
			    memcpy(seq[i], tempdummy, sizeof(char)*(lentotal));
			    fclose(fp);
			    return seq;
            }
			while(dummy[0]!='>'){
            if(dummy[0] == '\n' || dummy[0] == '\0'){
			    MALLOC(seq[i], char, lentotal+1);
			    memcpy(seq[i], tempdummy, sizeof(char)*(lentotal));
			    fclose(fp);
			    return seq;
            }                
			len = strlen(dummy);
			if(len == MAX_LINE) printf("There is very long line (> %u).\n", MAX_LINE);
			len += (dummy[len-1] == '\n')? -1:0;
			memcpy(tempdummy+lentotal, dummy, sizeof(char)*len);
			lentotal+=len;
            if(fgets(dummy, MAX_LINE, fp) == NULL){ /* read next fasta element */
			    MALLOC(seq[i], char, lentotal+1);
			    memcpy(seq[i], tempdummy, sizeof(char)*(lentotal));
			    fclose(fp);
			    return seq;
            }
            }
            MALLOC(seq[i], char, lentotal+1);
            memcpy(seq[i], tempdummy, sizeof(char)*(lentotal));
		}
		else if(dummy[0] == '\n' || dummy[0] == '\0') break;
		else printf("This is not a valid fasta file.\n");
	}
	if(fclose(fp)!=0){
         printf("Fail to close fasta file.\n");
    }
	return seq;
}

char** readBlockSeq(char *fname, unsigned int *num_seq)
{
	char dummy[MAX_LINE], **seq = NULL;
	unsigned int len = 0, lentotal = 0;
	*num_seq=0;
	fstream BlockSeq;
	BlockSeq.open(fname);
	if(BlockSeq.fail())
	    cout<<"Cannot open block sequence file"<<endl;
    while(!BlockSeq.eof())
    {
        BlockSeq.getline(dummy, MAX_LINE);
        if(dummy[0]=='>')
            (*num_seq)++;
    }
    if(*num_seq==0)
    {
        cout<<"Not A Valid Fasta File!"<<endl;
        return 0;
    }    
    BlockSeq.clear();   
    BlockSeq.seekg(0,ios::beg);
	cout<<"num_seq = "<<*num_seq<<endl;

    NewMem(*num_seq);
	MALLOC(seq, char*, *num_seq);
	MALLOC(SeqID, char*, *num_seq);
	EigenSeqs=new char**[*num_seq];

    for(int i=0;i<*num_seq; i++)
    {
        //**********Read Center Seq ID***************  
        BlockSeq.getline(dummy, MAX_LINE);
        len=strlen(dummy);
        len+=(dummy[len-1]=='\n')? -1:0;
        MALLOC(SeqID[i], char, len+1);
        memcpy(SeqID[i],dummy, sizeof(char)*len);
        //**************************************
        
        //**********Read Center Seq********************
        BlockSeq.getline(dummy, MAX_LINE);
        len=strlen(dummy);
        len+=(dummy[len-1]=='\n')? -1:0;
        MALLOC(seq[i], char, len+1);
        memcpy(seq[i],dummy, sizeof(char)*len);
        //**************************************
        
        BlockSeq>>Abundance[i];
        BlockSeq>>Sigma[i];
        BlockSeq.getline(dummy, MAX_LINE);
        int NumberOfSeqs=MIN(20,Abundance[i]);
        EigenSeqs[i]=new char*[NumberOfSeqs];
        
        //***********Read Eigen Seqs*************
	    for(int j=0;j<NumberOfSeqs;j++)
	    {
            BlockSeq.getline(dummy, MAX_LINE);
            len=strlen(dummy);
            len+=(dummy[len-1]=='\n')? -1:0;
            EigenSeqs[i][j]=new char[len+1];
            EigenSeqs[i][j][len]='\0';
            memcpy(EigenSeqs[i][j],dummy, sizeof(char)*len);
		}
		//*************************************************
	}
    BlockSeq.close();
	return seq;
} 

void NewMem(int nn)
{
     //Global Variables
     //*******************************
     Distrk=new double[Step];
     Time=new double[Step];
     Pi=new double[nn];
     PiPrev=new double[nn];
     Center=new int[nn];
     Sigma=new double[nn];
     SigmaChange=new bool[nn];
     Likelihood=new double[nn];
     LikelihoodRatio=new double[nn]; //for UpdateLikelihood;
     //************************************
     LikelihoodTemp=new double*[nn];
     for(int i=0;i<nn;i++)
         LikelihoodTemp[i]=new double[1000];
     //************************************
     X=new float*[nn];         //Distance Matrix
     for(int i=0;i<nn;i++)
         X[i]=new float[nn];
     Z=new int[nn];       //missing data
     NN=new int[nn];       //NN[i-1]=# of samples in i-th cluster
     IterNN=new int[nn];
     Abundance=new int[nn];
     Centerlikelihood=new double[nn];  //square sum of current center
     Weight=new bool[nn];  //for Birth-Death control;
     //*******************************

     //Result Process
     //*******************************
     ResultPi=new double[nn];
     ResultCenter=new int [nn];
     ResultSigma=new double[nn];
     ResultLikelihood=new double[nn];
    //*******************************

    //Hyperparameters
    //*******************************
    //For Pi
    Gamma=new double[nn];   
    return;
}

void DeleteMem(int nn, int iter)
{
     //Global Variables
     //*******************************
     DELETE(Distrk);
     DELETE(Time);
     DELETE(Pi);
     DELETE(PiPrev);
     DELETE(Center);
     DELETE(Sigma);
     DELETE(SigmaChange);
     DELETE(Likelihood);
     DELETE(LikelihoodRatio);
     for(int i=0;i<nn;i++)
     {
             FREE(SeqID[i]);
             FREE(seq[i]);
     }
     FREE(SeqID);
     FREE(seq);
     if(iter!=0)
     {
     for(int i=0;i<nn;i++)
     {
        for(int j=0;j<MIN(20,Abundance[i]);j++)
            delete[] EigenSeqs[i][j];
        delete[] EigenSeqs[i];
     }
     delete[] EigenSeqs;
     }   
     //for UpdateLikelihood;
     //************************************
     for(int i=0;i<nn;i++) 
             delete[] LikelihoodTemp[i];
     delete[] LikelihoodTemp;
     //************************************
     for(int i=0;i<nn;i++)
	         delete[] X[i];
     delete[] X;         //Distance Matrix
     DELETE(Z);       //missing data
     DELETE(NN);       //NN[i-1]=# of samples in i-th cluster
     DELETE(IterNN);
     DELETE(Abundance);
     DELETE(Centerlikelihood);  //square sum of current center
     DELETE(Weight);  
     //for Birth-Death control;
     //*******************************

     //Result Process
     //*******************************
     DELETE(ResultPi);
     DELETE(ResultCenter);
     DELETE(ResultSigma);
     DELETE(ResultLikelihood);
    //*******************************

    //Hyperparameters
    //*******************************
    //For Pi
    DELETE(Gamma);
    return;
}
